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abstract 


This  report  concerns  the  use  of  frequency  domain  computer  codes  such  as  the 
Numerical  Electromagnetic  Code  (NEC)  for  computing  the  time  domain  Electro¬ 
magnetic  Pulse  (EMP)  response  of  structures  such  as  antennas,  aircraft  or 
communication  shelters.  The  proper  representation  of  the  EMP  excitation  and  the 
selection  of  a  number  of  appropriate  frequencies  to  obtain  a  correct  time  domain 
EMP  response  are  studied.  The  effects  of  adapting  the  modelling  of  the  problem 
for  different  frequency  ranges  Is  discussed.  Guidelines  are  given  for  obtaining 
a  correct  time  domain  response  with  efficient  use  of  computer  time. 


Ce  rapport  consldAre  1 'utilisation  de  programmes  calculant  la  reponse  en 
frequence  de  syst^mes  tels  que  des  antennes,  des  avions  ou  des  systemes  de 
communications  pour  obtenlr  la  rdponsa  temporelle  resultant  dune  excitation  de 
type  lEM  (Impulsion  dlectromagndtlque) .  Les  consequences  dune  representation 
approprie  de  1  excitation  lEM  alnsl  que  de  la  selection  dun  certain  nombre  de 
frequences  pour  calculer  la  transformee  de  Fourier  inverse  sont  etudiees. 
Quelques  regies  simples  sont  donndes  pour  obtenir  une  reponse  temporelle  exacte 
tout  en  utllisant  le  minimum  de  temps  calcul. 


ill 


The  Electromagnetic  Pulse  (EMP)  interaction  analysis  Is  an  important  part 
of  the  recommended  procedure  for  hardening  electronic  equipment  against  the 
effects  of  EMP  generated  by  a  nuclear  detonation.  This  report  concerns  the  use 
of  frequency  domain  computer  codes  such  as  the  Numerical  Electromagnetic  Code 
(NEC)  for  computing  the  time  domain  EMP  response  of  structures  such  as  antennas, 
aircraft  or  communication  shelters. 

Some  properties  of  the  Discrete  Fourier  transform  (DFT)  and  the  Fast 
Fourier  transform  (FFT)  and  some  of  their  slde>effects  when  applied  to  solve  time 
domain  problems  are  discussed.  A  technique  using  cublc>splina  interpolation  to 
perform  a  DFT  or  a  FFT  on  unevenly  spaced  sequences  is  introduced. 

The  effects  of  the  choice  of  a  proper  model  for  the  EMP  excitation  and  of 
the  selection  of  an  appropriate  number  of  frequencies  for  taking  an  inverse 
Fourier  transform  for  obtaining  a  correct  time  domain  EMP  response  are  studied. 
Comparisons  are  made  between  the  results  obtained  by  the  Fourier  transformation 
of  the  frequency  domain  results  and  chose  obtained  with  a  time  domain  code  such 
as  the  Thin  Wire  Time  Domain  code  (TWTD) .  General  guidelines  are  given  for 
obtaining  a  correct  time  domain  response  with  efficient  use  of  computer  time. 
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1.0  INTRODOCTION 


This  report  concerns  the  use  of  frequency  domain  computer  codes  such  as  the 
Numerical  Elec tromag'ne tic  Code  (NEC)  for  computing  the  time  domain  Electro¬ 
magnetic  Pulse  (EMP)  response  of  structures  such  as  antennas,  aircraft  or 
communication  shelters.  The  effects  of  the  choice  of  a  proper  model  for  the 
excitation  and  of  the  selection  of  an  appropriate  number  of  frequencies  for 
taking  an  inverse  Fourier  transform  and  obtaining  a  correct  time  domain  EMP 
response  are  studied.  Guidelines  are  given  for  obtaining  a  correct  time  domain 
response  with  efficient  use  of  computer  time.  Validation  of  this  technique  is 
obtained  by  comparing  the  results  of  a  frequency  domain  code  with  chose  of  a  time 
domain  code  as  well  as  with  experimental  results. 


1.1  EMP  ANALYSIS  AND  HARDENING 

Electromagnetic  pulse  resulting  from  a  nuclear  detonation,  usually  referred 
to  as  EMP,  is  one  of  Che  most  sericus  nuclear  effects,  since  it  can  adversely 
affect  performances  of  electronic  and  weapon  systems  hundreds  of  kilometres  away 
from  Che  source.  The  EMP  is  a  high -intensity,  short -duration  electromagnetic 
field  and  although  its  energy  content  is  not  very  large  because  of  its  short 
duration,  electronic  components  or  systems  can  be  upset  or  permanently  damaged 
by  the  very  high  voltages  or  currents  which  can  be  induced  and  coupled  into  the 
system.  Analytical  cools  and  methodology  have  been  developed  to  assess  EMP 
vulnerability  and  obtain  EMP  hardness  of  a  system.  Different  EMP  simulators  and 
test  procedures  have  also  been  developed  to  ensure  EMP  hardness. 

Hardening  against  EMP  is  an  iterative  process  where  the  EMP  interaction 
with  Che  system  is  predicted,  the  effectiveness  of  the  shielding,  filters,  and 
other  protective  components  is  estimated,  and  induced  voltages  or  currents 
reaching  sensitive  components  are  checked  against  their  upset  or  damage 
threshold.  EMP  testing  is  performed  to  increase  confidence  chat  system  hardening 
has  been  properly  addressed.  Susceptibilities  discovered  during  analysis  or 
testing  are  corrected  with  hardening  techniques,  such  as  shielding  and  filtering. 


1.2  EMP  NUMERICAL  ANALYSIS  CODES 

The  complete  problem  of  Che  coupling  of  the  EMP  frcm  a  specific  detonation 
into  Che  components  of  a  specific  system  is  an  extremely  complex  problem.  There 
are  no  existing  computer  codes  which  can  treat  the  complete  problem,  and  as  a 
consequence,  various  aspects  of  Che  EMP  problem  are  studied  individuallv  with  the 
aid  of  Che  existing  codes.  Good  engineering  practices  must  then  be  used  to 
combine  the  various  parts  of  Che  solution  into  an  assessment  of  the  E.MP 
vulnerability. 

The  various  computer  codes  available  can  be  divided  into  three  general 
classes:  environment  codes,  coupling  codes,  and  circult-analvsls  codes.  The 
environment  codes  calculate  the  electromagnetic  fields  generated  by  a  nuclear 
detonation.  The  coupling  codes  determine  the  currents,  potentials,  and  charges 
induced  into  isolated  bodies  and  transmission  lines  by  incident  electromagnetic 
fields.  The  circuit-analysis  codes  analyze  Che  response  of  linear  and  nonlinear 
circuits  to  voltage  and  current  injections  resulting  from  EMP  illumination. 
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Tha  anvlronaant  codas  ara  usad  Co  pradlec  tha  alaccromagnecic  fields 
ganaracad  by  nuclaar  daconaclons  and  visually  raqulra  in>dapch  knowledge  of  a 
weapon.  For  practical  anginearlng  practices,  tha  EMF  environment  is  usually 
specified  for  various  type  of  equipments  or  installations  in  standards  documents, 
such  as  [ 2 ] . 

The  alactromagnetlc  interaction  coupling  codas  are  used  to  calculate  the 
charges  and  currents  Induced  in  structures  or  transmission  lines  by  an  Incident 
pulse  and  estimate  tha  field  penetration  through  apertures  or  shielding. 
Designers  of  those  codes  have  taken  different  mathematical  or  numerical 
approaches  [S]  which  make  them  better  suited  for  certain  type  of  structures  or 
problems.  Tha  interaction  codes  can  be  divided  into  two  general  categories: 
^  frequency  domain  codas  and  time  domain  codes.  Frequency  domain  codes  evaluate 

yr,  the  currents  and  charges  induced  at  a  specified  frequency.  A  large  number  of 

/  frequencies  is  required  in  order  to  estimate  the  time  domain  response  of  the 

system  from  its  frequency  respon.se.  Tima  domain  codes  compute  the  time  response 
directly,  which  make  them  more  suitable  to  solve  EHP  problems.  However,  in  many 
cases,  a  time  domain  code  suitable  to  solve  a  particular  problem  does  not  exist 
and  a  frequency  domain  code  must  be  usad. 

Tha  cireuit<analyals  programs,  such  as  SPICE  or  MICRO*CAP,  are  primarily 
nonlinear,  translent'analysis  codes.  Some  will  also  perform  AC.  DC,  and 
frequency* analysis  calculations.  They  do  not  compute  the  coupling  of  the  EMP 
transient  fields  into  the  circuit  directly  and  thus  rely  on  other  codes  to 
provide  the  EMP  coupling,  usually  as  voltages  or  currents  as  functions  of  time. 
Performances  of  protective  devices,  such  as  arrestors  and  filters,  can  be 
evaluated  to  predict  tha  currents  or  voltages  reaching  sensitive  components  of 
a  system. 


1.3  ELECTROMAGNETIC  INTERACTION  CODES 

This  section  briefly  describes  several  interaction  codes  chat  are  currently 
in  use  at  DREO  or  being  implemented.  A  frequency  domain,  NEC,  and  a  time  domain 
code,  TUTD,  were  chosen  to  conduct  this  study  to  obtain  comparative  data  between 
Che  two  methods.  All  computer  simulations  were  executed  on  a  VAX  11/780  or  on 
a  MicroVAX-II.  All  execution  times  (CPU  time)  given  are  normalized  for  a 
MlcroVAX-II. 


1.3.1  Thin  Wire  Tima  Domain  Code 

The  Thln-Wlre  Time-Domain  coda  (TWTD)  [6]  is  a  time  domain  computer  code 
CO  compute  Che  induced  currents  on,  and  the  radiated  or  scattered  fields  from  an 
arbitrary  thin  wire  antenna  or  structure.  The  excitation  on  the  structure  is  a 
specified  time-varying  source  or  an  incident  plane  wave.  The  output  may  include 
currents,  the  radiated  or  scattered  fields,  antenna  gain,  and  the  spectral 
characteristics  of  Che  input  impedance.  Although  TWTD  analyzes  only  structures 
in  free  space,  structures  over  perfect  ground  can  be  analyzed  by  explicitly  using 
image  theory.  The  code  was  modified  to  accommodate  larger  structures  or  longer 
time  history  and  to  define  new  excitations  specific  to  EMP  studies. 


2 


1.3.2  Numerical  Electromagnetic  Code 


The  Numerical  Electromagnetic  Code  (NEC)  [7]  is  a  frequency  domain  computer 
code  for  the  analysis  of  the  electromagnetic  response  cf  antennas  and  other  metal 
structures.  The  code  combines  formulation  of  smooth  surfaces  and  wires  for 
convenient  and  accurate  modelling  of  a  wide  variety  of  structures.  A  model  may 
Include  nonradiating  networks  and  transmission  lines  connecting  parts  of  the 
structure,  perfect  or  imperfect  conductors,  and  lumped-element  loading.  A 
structure  may  also  be  modeled  over  a  ground  plane  that  may  be  either  a  perfect 
or  a  imperfect  conductor.  The  excitation  may  be  either  voltage  sources  on  the 
structure  or  an  incident  plane  wave.  The  output  may  include  the  induced  currents 
and  charges,  the  near  electric  or  magnetic  fields,  and  the  radiated  fields. 

NEC  is  probably  one  of  the  most  versatile  code  available  today.  NEC  has 
been  used  for  many  radiation  and  scattering  problems.  Its  use  for  EMP  coupling 
problems  has  been  less  common  because  of  the  large  amount  of  computing  time 
required  for  the  determination  of  currents  on  a  complex  structure  for  a  large 
number  of  frequencies. 


1.3.3  Electric-Field  Integral  Equation  Code 

The  Electric-Field  Integral  Equation  Code  (EFIE)  [8]  is  a  frequency  domain 
computer  code  developed  to  treat  the  problem  of  elactromagnetic  scattering  from 
arbitrarily  shaped,  perfectly  conducting  objects,  illuminated  by  an  incident 
plane  wave.  Objects  are  modelled  with  triangular  surface  patches  which  may 
include  apertures  or  cavities.  The  program  computes  the  surface  current  density 
induced  on  either  an  open  or  closed  surface  as  well  as  the  bistatic  radar  cross 
section  of  the  object.  The  code  was  modified  in  house  to  compute  also  the  field 
intensities  from  the  current  densities  on  the  surface.  There  have  been  other 
major  modifications  made  to  produce  drastic  reduction  in  computing  time. 


1.3.4  General  Electromagnetic  Model  for  Analysis  of  Complex  Systems 

The  General  Electromagnetic  Model  for  Analysis  of  Complex  Systems  (GEM<\CS) 
incorporates  a  variety  of  techniques  for  electromagnetic  analysis  of  complex 
objects.  The  latest  available  version  of  the  code  includes  thin  wire  and  surface 
patch  Method  of  Moments  (MOM)  formalism,  with  or  without  Geometrical  Theory  of 
Diffraction  (GTD)  interactions,  to  solve  for  exterior  problems.  It  also  includes 
a  Finite  Difference  (FD)  formalism  in  the  frequency  domain  to  solve  interior 
problems  and  mathematics  necessary  to  connect  exterior  and  interior  solutions 
when  apertures  are  present. 
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2.0  IHOISCREIB  J9PRIE&..TRAMSr9RM 


Th«  procedure  described  in  chls  report  is  used  to  solve  time  domain 
problems.  It  is  based  essentially  on  computing  the  frequency  domain 
representation  of  the  solution  and  applying  an  inverse  discrete  Fourier  transform 
algorithm  to  obtain  the  time  domain  solution.  This  chapter  introduces  the 
general  procedure  used  to  obtain  the  solution  and  provides  a  detailed  description 
of  discrete  Fourier  transform  and  its  properties. 


2.1  GENERAL  PROCEDURE 

A  time  domain  analysis  code,  if  one  is  available  to  solve  a  particular 
problem,  is  usually  more  efficient  to  compute  time  domain  solutions.  In  many 
cases  however,  the  only  choice  is  to  use  a  frequency  domain  analysis  code.  Many 
frequency  domain  codes  have  been  developed  to  solve  radiation  and  scattering 
problems.  Some  of  them,  such  as  NEC,  are  very  versatile  and  can  be  used  to  model 
a  large  variety  of  problems. 

A  linear  time* invariant  system  can  be  fully  characterized  by  its  transfer 
function  H((i)) .  A  frequency  domain  analysis  code  can  be  run  at  multiple 
frequencies  to  compute  H((j) .  The  response  y(t)  of  a  system  to  a  known  excitation 
e(c)  can  be  expressed  in  frequency  domain  as: 

Y(w)  -  H(w)  .  E<«)  (1) 


where  ECw)  and  Y(w)  are  the  Fourier  transform  of  e(c)  and  y(c)  respectively  and 
H(w)  is  the  numerically  computed  transfer  function.  For  our  purpose,  the 
excitation  e(t)  is  known  and  its  Fourier  transform  can  be  evaluated  analytically. 
The  solution  y(c)  is  obtained  with  the  inverse  Fourier  transform  of  H(w)E(w): 

e(t)  •  E(w) 

H(u)  (2) 

U 

y{t)  •  Y(w) 

This  procedure  Involves  using  a  frequency  domain  code  such  as  NEC  to 
compute  the  frequency  response  H(w) .  This  is  very  computationally  intensive, 
easily  taking  several  days  of  CPU  time  to  complete.  In  addition,  the  transfer 
function  is  not  computed  in  its  continuous  form,  but  rather  at  a  number  of 
discrete  frequencies.  Consequently,  the  discrete  form  of  the  Fourier  transform 
(DFT)  must  be  used.  Some  of  the  DFT  properties  must  be  understood  to  avoid  some 
'side-effects'  when  used.  Guidelines  to  minimize  the  number  of  frequencies  and 
CPU  time  required  will  be  presented  in  Chapter  3 . 


2.2  THE  DISCRETE  FOURIER  TRANSFORM  ALGORITHM 

The  Fourier  transform  (FT)  of  a  continuous  waveform  is  a  continuous - 
frequency  signal  representation  of  its  spectrum  or  frequency  content.  However, 
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when  a  waveform  Is  sampled,  which  is  necessary  when  Che  system  is  Co  be  analyzed 
on  a  computer,  its  discrete  counterpart,  the  discrete  Fourier  transform  (DFT) 
must  be  used.  The  DFT  retains  most  of  Che  properties  of  Che  continuous  Fourier 
transforms.  There  are  however  some  important  differences  which  may  cause  the  DFT 
to  produce  erroneous  results.  Some  properties  more  relevant  to  our  application 
will  be  discussed  here.  Many  text  books  treat  Che  FT  and  Che  DFT  in  more  detail 
(16]  [11]  [12]  [13]. 

The  Fourier  transform  and  its  inverse  are  defined  as: 


X<f)  -  f  x(c)  e--12*«dt  (3) 


and 


x(c) 


£  X(f)  eJ^'fftdf 


W 


As  observed  in  equations  (3)  and  (4),  some  symmetry  or  duality  exists 
between  the  two  equations.  This  duality  can  be  formulated  as: 

if  x(t)  -  X(f) 

then  X(c)  ••  x(-f) 

This  duality  implies  that  every  property  of  Che  Fourier  transform  is  also 
applicable  to  its  inverse.  It  follows  that  if  there  are  characteristics  in  time 
domain  which  have  implications  in  frequency  domain,  then  Che  same  characteristics 
in  frequency  domain  will  have  similar  implications  in  time  domain. 

The  discrete  form  of  the  Fourier  transform  (DFT)  and  its  inverse  (IDFT)  are 
defined  as^: 


N-l 

X(k)  -  x(n)e--'2’'>'n/N  (6) 

n«0 

and 

x(n)  -  4  y  X(k)ej2’"‘'’'''  (7) 

where  x(n)  is  a  real  sequence^  of  N  points  representing  a  waveform  sampled  at 
a  regular  interval  At  and  X(k)  a  complex  sequence  of  N  points  representing  its 
transform  at  frequencies  0,  fo,  2fo.  ...  ,  (N-Dfg.  It  is  clear  from  these 
definitions  that  both  sequences  X(n)  and  X(k)  are  periodic  with  a  period  of  N 
samples.  A  consequence  of  this  periodicity  is  that  shifting  a  sequence  always 
implies  a  circular  shift.  The  effect  of  this  property  in  time  domain  is  chat 


^  The  definition  of  the  DFT  is  not  uniform  in  the  literature.  The  1/N 
factor  is  sometime  placed  in  the  X(k)  definition. 

^  In  general,  the  sequence  x(n)  can  be  complex.  In  our  case  however,  the 
solution  is  known  to  be  real  as  it  represent  a  measurable  quantity 
(current  or  voltage). 
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whsn  shifting  a  sequenca,  the  end  of  the  wavefom  wraps-around  Its  beginning. 
The  effect  in  frequency  domain  Is  chat  the  sequence  appears  to  have  been  mirrored 
around  N/2  when  in  fact  the  frequencies  at  k-N  to  N/2  are  the  negative 
frequencies  of  the  transform  from  >df  Co  » f,../2 . 

It  can  be  derived  from  the  symmetry'  property  of  Che  Fourier  transform  chat 
if  x(n}  is  real,  then 

X(k)-X*(-k)  (8) 


where  X*  denotes  Che  conjugate.  It  states  that  before  a  time  domain  sequence  can 
be  reconstructed,  the  frequency  response  obtained  from  the  simulations  must  be 
first  mirrored  with  its  conjugate. 

The  sampling  interval  At  and  Che  fundamental  frequency  fo  or  Af  do  not 
appear  explicitly  in  aquations  (6)  and  (7).  The  time  and  frequency  scaling 
property  of  the  Fourier  transform 

x(at)-^G(i)  (9) 


can  be  applied  to  adjust  the 
sequences.  It  follows  Chat; 

amplitude  and 

calculate  the 

interval  of  Che 

Ac  •  -  and  Af  •  .,^11 . 

TTaT  TTSt 

(10) 

and 

t^  -  (N-l)At  «  ^  and 

-  (N-l)Af  » 

1 

AC 

(11) 

For  EHP  problems,  a  finite  duration  excitation  yields  a  finite  duration 
response.  Consequently,  Che  frequency  sampling  Af  must  be  small  enough  to 
translate  into  a  sufficient  time  duration  t,.,.  Similarly,  At  must  be  small 
enough  to  capture  the  fastest  transitions  of  Che  solution,  imposing  a  sufficient 
frequency  range  f„„,  or  for  a  given  Af,  a  sufficient  number  of  frequencies  N. 

A  common  technique  used  with  the  DFT  consists  in  extending  a  sequence  x(n) 
of  N  points  Co  M  points  by  adding  zero-valued  samples  at  the  end  (zero-filling). 
It  can  be  easily  demonstrated  from  equations  (6)  and  (7)  that  this  allows  one  to 
compute  Che  transform  at  arbitrary  resolution,  effectively  yielding  exact 
interpolation  between  points.  This  property  is  most  useful  to  fit  Che  length  of 
Che  sequence  to  meet  the  limitations  of  some  of  the  practical  DFT 
implementations . 


2.3  THE  FAST  FOURIER  TRANSFORM  ALGORITHM 

The  fast  Fourier  transform  (FFT)  is  merely  a  very  clever  implementation  of 
Che  DFT.  For  a  sequence  of  N  points,  the  execution  time  is  proportional  to  4N^ 
for  the  DFT  compared  to  N/2*log2N  for  Che  FFT.  For  a  1024  points  sequence,  the 
FFT  algorithm  is  over  800  times  faster.  The  original  FFT  implementation  [14]  was 
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chosen,  with  few  modifications  to  take  into  account  the  scaling  due  to  df  as 
defined  by  equation  (9).  A  listing  of  the  subroutine  Is  given  In  Appendix  A. 
The  only  restriction  with  this  algorithm  Is  that  the  size  of  the  sequence  Is 
restricted  to  powers  of  2  (N-2") .  A  more  exhaustive  overview  of  alternate 
Implementation  of  the  FFT  can  be  found  In  [15] . 


2.4  RELATIONSHIPS  BETWEEN  THE  DFT  AND  THE  FOURIER  TRANSFORM 

Whenever  the  lOFT  Is  used  to  reconstruct  a  continuous  waveform  from  its 
frequency  domain  definition,  a  set  of  operations,  whether  they  are  Implied  or 
explicit,  are  performed,  each  altering  the  result.  To  illustrate  the  most  common 
problems  with  the  DFT  (aliasing,  leakage,  etc.),  consider  a  waveform  y(t)  derived 
from  its  known  Fourier  transform  Y(f),  as  shown  on  Figure  2>l(a). 

The  transform  Y(f)  Is  computed  from  ntimerlcal  simulations,  at  a  number  of 
equally  spaced  frequencies.  This  corresponds  to  multiplying  Y(f)  with  a  Dirac 
comb  function,  with  a  period  of  Af.  The  transform  to  time  domain  of  this  comb 
function  is  also  a  comb  function  with  a  period  of  1/Af  (Figure  2-l(b)). 
Multiplication  in  frequency  domain  yields  a  convolution  in  time  domain,  resulting 
in  y(t)  becoming  periodic  and  being  corrupted  with  aliasing  (Figure  2*l(c)) .  The 
only  method  to  reduce  the  aliasing  Is  to  decrease  Af . 

Only  a  finite  portion  of  Y(f)  Is  computed,  for  frequencies  from  0  to  f„,j(, 
thus  ignoring  some  higher  frequencies.  This  is  equivalent  to  multiplying  Y(f) 
with  a  rectangular  window  W(f),  of  width  ±f„...  As  observed  on  Figure  2-l(d), 
Che  rectangular  window  has  a  transform  in  the  form  of  a  ain(x)/x  function.  When 
Y(f)  is  convolved  with  the  window's  transform,  some  leakage  effect  can  be 
observed  (Figure  2-l(a)).  This  effect  can  be  reduced  by  increasing  f„g,,  and 
thus  N,  or  the  using  a  different  window,  such  as  a  cosine- capered  window.  The 
combined  effect  of  the  sampling  and  windowing  is  shown  on  Figure  2-l(f). 
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3.0  SOLVING  TIME  DOMAIN  PROBLEMS  WITH  THE  DFT 


The  procedure  described  Co  obtain  a  time  domain  solution  by  caking  the 
inverse  discrete  Fourier  transform  of  a  computed  frequency  domain  solution  is 
relatively  simple.  However,  in  Che  absence  of  any  guidelines,  this  process  can 
be  very  time  consuming  and  lead  to  erroneous  results.  Simple  guidelines  will  be 
developed  in  this  chapter  to  obtain  an  accurate  time  domain  solution  with  the 
efficient  use  of  computer  time. 

To  illustrate  some  of  Che  techniques  presented  in  this  chapter,  a  sample 
problem  will  be  used.  Figure  3*1  shows  a  stick  model  of  a  helicopter.  It  is  a 
crude  wire  grid  representation  of  it,  consisting  of  185  segments.  The  currents 
induced  from  an  incident  EMP  waveform  were  computed  with  a  time  domain  code 
(TMTD)  and  with  a  frequency  domain  code  (NEC).  In  this  example,  it  took  NEC  over 
136  hours  on  a  MicroVAX  II  to  compute  Che  currents  at  1024  frequencies,  0.25  MHz 
apart,  giving  the  frequency  response  of  the  system  as  shown  on  Figure  3-2, 
necessary  to  obtain  an  accurate  time  domain  response,  shown  on  Figure  3-3.  It 
Cook  TUTD  about  2.5  hours  to  compute  a  solution,  but  due  to  code  limitation,  only 
the  first  400  ns  were  obtained.  Nevertheless,  the  agreement  with  NEC  is  very 
good,  as  shown  on  Figure  3-4.  This  solution  will  be  used  thorough  this  chapter 
as  a  point  of  comparison  to  evaluate  Che  effect  of  various  parameters. 


3.1  EMP  WAVEFORM  DEFINITION 

The  EMP  produced  by  a  nuclear  burst  at  high  altitude  is  a  large-amplitude, 
very  short  duration  transient  field  covering  a  very  wide  area  beneath  the  burst 
point.  In  general,  the  exact  characteristics  of  the  EMP  field  such  as  peak 
amplitude,  rise  time  and  polarization,  depend  on  many  factors,  such  as  Che  weapon 
yield,  the  height  of  burst  and  the  observer's  location.  For  practical  purposes, 
a  standard  waveform  representing  a  worst  case  is  used  [2]  [3]. 

The  simplest  analytical  expression  approximating  the  EMP  waveform,  often 
referred  to  as  the  Bell  curve  [4]  or  the  old  NATO  definition,  is  the  double 
exponential  waveform; 

E(C)  -  AV  .(e-<*‘-e-®‘)  (12) 


As  shown  in  Figure  3-5,  this  pulse  has  a  peak  value  of  50  kV  per  meter,  a 
rise  time  (10-90%)  of  about  4.2  nsec,  a  pulse  width  (50-50%)  of  about  185  nsec 
and  a  decay  time  (peak-to-10%)  of  about  600  nsec.  This  expression  is  very  useful 
due  to  its  simplicity.  Its  Fourier  transform: 


E(w)-AV  ■  (  -} . . 

a+«j 


(13) 


and  its  Laplace  transform  are  very  simple,  which  makes  it  possible  to  solve 
simple  problems  analytically.  However,  this  expression  has  a  major  flaw  which 
can  render  the  results  of  a  numerical  simulation  useless.  As  observed  in 
Figure  3-5,  the  waveform  is  continuous  over  time,  but  its  first  derivative  is 
discontinuous  at  t-t(j.  Since  the  coupling  of  electromagnetic  waves  to  a 
structure  is  strongly  dependent  of  Che  derivative  of  the  excitation,  this 
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Figure  3-1  Stick  model  of  a  helicopter,  consisting  of  185  thin  wire  segments. 

The  model  is  about  19.5  meters  long.  Calculations  were  made  with 
NEC  and  TWTD  to  estimate  the  current  at  the  nose  (segment  on  the 
right) . 
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Figure  3-2  Transfer  function  (magnitude  and  phase)  at  the  nose  of  the  stick 
model  helicopter,  computed  by  NEC  at  1024  frequencies,  0.25  MHz 


CURRENT  (A) 


CURRENT  (A) 


Figure 


TIME  (ns) 


3-4  Current  at  the  nose  of  the  stick  model  helicopter  induced  from  an 
incident  EMP  waveform  calculated  with  NEC  (slightly  shifted) 
compared  with  the  current  calculated  with  TWTD. 


discontinuity  can  introduce  significant  ringing  and  distortion  in  the  solution. 

A  better  expression,  known  as  the  new  NATO  definition,  to  approximate  the 
EHP  waveform  is: 


E(t) 


AU  m 


(U) 


As  shown  in  Tigure  3-3,  it  has  the  same  characteristics  as  the  Bell  curve 
(peak  amplitude,  rise  time  and  pulse  width),  but  its  derivatives  are  now 
continuous  over  time  and  it  has  a  more  realistic  leading  edge.  Note  that  the  AV, 
a  and  0  parameters  are  different  from  (12).  The  Fourier  transform  of  (14)  is: 


E(u) 


AV 


It 


sin()r* 


) 


(15) 


Figure  3-6  shows  the  frequency  spectrum  for  both  definitions  as  calculated 
with  (13)  and  (IS).  The  most  Important  feature  is  that  most  of  the  energy  is 
concentrated  below  100  KHz  with  virtually  nothing  left  above  500  KHz. 


3.2  MINIMIZING  THE  COMPUTATION  TIME  FOR  FREQUENCY  DOMAIN  CODES 

In  order  to  minimize  the  CPU  time,  two  steps  can  be  taken:  reduce  the 
number  of  segments  or  reduce  the  nvuober  of  frequencies.  For  instance,  NEC 
execution  time  grows  exponentially  with  the  number  of  segments  or  patches  and 
consequently,  keeping  this  number  as  low  as  possible  will  translate  in 
considerable  savings  in  CPU  time.  NEC  guidelines  state  that  segment  length 
should  be  smaller  chan  1/10  at  Che  highest  frequency.  By  adjusting  the  segment 
length  for  the  different  frequencies,  i.e.  use  longer  segments  at  lower 
frequencies,  it  is  possible  to  save  up  to  50%  of  CPU  time.  For  example, 
modelling  the  helicopter  by  using  a  63  segment  model  for  frequencies  up  to  64  MHz 
and  a  185  segment  model  for  frequencies  up  to  256  MHz  gives  a  saving  of  22%. 
However,  when  varying  Che  number  of  segments,  it  is  important  to  keep  Che  point 
of  observation  constant  and  to  carefully  look  at  the  data  to  ensure  continuity 
when  passing  from  one  model  to  another. 


3.3  SELECTION  OF  AN  ADEQUATE  FREQUENCY  STEP 

The  two  important  parameters  chat  determine  the  number  of  frequencies  are 
the  frequency  step  (Af)  and  frequency  range  (fmax)-  Figure  3-7  shows  the  effect 
of  various  frequency  steps  used  in  obtaining  the  time  domain  response. 
Obviously,  if  Af  is  chosen  coo  large,  some  Important  features  of  the  frequency 
response  such  as  Che  resonance  peaks  may  be  missed.  As  seen  in  equation  (11), 
there  is  a  direct  relation  between  the  frequency  step  and  the  time  duration  of 
Che  solution  inherent  to  the  inverse  Fourier  transform.  If  the  Af  is  chosen  too 
large,  the  time  duration  of  the  solution  (t„„)  will  be  shorter  chan  the  response 
duration  and  Che  effect  will  be  chat  the  solution  will  not  decay  to  zero.  This 
can  be  observed  in  Figure  3-7  where  the  4  MHz/step  curve  has  a  250  ns  duration 
and  is  obviously  incomplete.  In  our  example,  the  response  duration  is  about  2  ms 


16 
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Figure  3-6  Frequency  spectrum  of  old  (double -exponential)  and  new  (reciprocal) 
NATO  definitions  for  the  EMP  waveform. 


and  therefore,  a  0.5  MHz  frequency  step  ts  adequate.  It  Is  also  worth  noting 
that  there  is  no  significant  improvement  when  going  from  0.3  MHz  to  0.23  MHz 
step.  An  initial  estimate  of  Af  can  be  obtained  by  estimating  the  first 
resonance  of  the  structure  from  its  the  shape  and  dimensions.  Starting  with  Af 
approximately  equal  to  1/10*’‘‘  of  the  estimated  first  resonance  frequency  is 
usually  adequate.  The  length  of  the  helicopter  for  instance  is  19.3  m,  which 
corresponds  to  a  resonance  around  8  MHz.  A  first  estimate  of  the  solution  can 
be  made  by  calculating  the  currents  at  1  MHz  steps. 


3.4  SELECTION  OF  AN  ADEQUATE  FREQUENCY  RANGE 

The  frequency  domain  simulation  codes  compute  the  solution  at  N 
frequencies,  from  Af  to  N-Af^,  effectively  ignoring  the  higher  harmonics  of  the 
solution.  A  reduced  frequency  range  results  in  more  distortion  (slower  rise 
time,  overshoot,  etc.)  due  to  the  fewer  number  of  harmonics  included  in  the 
solution.  It  is  easier  to  estimate  what  the  proper  range  should  be  by  plotting 
the  solution  H(w)'E(w)  in  the  frequency  domain  as  shown  of  Figure  3-8.  From  even 
a  limited  number  of  frequencies,  calculated  at  a  possibly  larger  Af,  Ic  is 
possible  to  observe  a  tendency  in  the  curve  indicating  chat  the  higher 
frequencies  contribute  less  and  less  to  the  solution.  This  is  especially  true 
for  frequencies  above  100  MHz  where  the  frequency  content  of  the  excitation  drops 
dramatically.  Figure  3-9  shows  Che  effect  of  the  frequency  range  on  the 
solution.  A  range  of  32  MHz  is  obtained  by  rejecting  the  higher  frequencies 
which  contribute  less  chan  10%  of  the  peak  resonance.  This  may  be  insufficient 
if  a  good  estimation  of  Che  rise  time  or  any  sudden  changes  is  desired,  but  may 
still  provide  an  approximation  of  the  solution.  If  the  range  is  extended  to  keep 
Che  frequencies  which  contribute  more  chan  1%  of  the  peak,  or  about  100  MHz,  Che 
solution  becomes  much  more  accurate  and  only  a  slight  overshoot  and  a  slight 
distortion  is  observed. 


3.5  INTERPOLATION  OF  THE  FREQUENCY  RESPONSE 

It  was  shown  that  Che  selection  of  an  appropriate  frequency  step  and 
frequency  range  has  a  major  impact  on  the  accuracy  of  the  response.  Ic  is  also 
evident  from  the  frequency  domain  response  (Figure  3-2)  chat  a  smaller  step 
should  be  used  around  Che  resonance  peaks,  especially  around  the  first  few  peaks. 
However  Che  standard  application  of  the  FFT  algorithm  requires  chat  all 
frequencies  be  equally  spaced,  thus  we  would  need  to  compute  the  currents  at  the 
same  spacing  for  the  entire  frequency  range  resulting  in  an  increase  of  the  CPU 
time  by  Che  same  factor. 

To  avoid  this  excessive  use  of  CPU,  a  code  was  developed  which  Cakes  a 
spectrum  H((j)  consisting  of  unevenly  spaced  frequencies  and,  by  using  cubic 
spline  interpolation,  generates  a  new  spectrum  at  evenly  spaced  frequencies  at 
a  different  and  possibly  a  smaller  frequency  step.  This  code  allows  us  to  use 
a  smaller  step  only  near  the  resonance  peaks  where  the  magnitude  and/or  phase 
change  rapidly.  Calculations  in  ocher  regions  can  be  made  at  a  bigger  step. 


^  The  solution. at  DC  is  always  0.  as  the  response  to  an  incident  field  is 
function  of  its  derivative. 
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Figure  3-7  Effect  of  the  frequency  step  on  the  time  domain  computed  solution 
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Figure  3-8  Produce  of  the  computed  transfer  function  (H(w))  and  the  spectrum  of 
the  EMP  excitation  (E(w)). 
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Figure  3-9  Effect  of  the  frequency  range  on  the  time  domain  computed  solution 
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Spline  incerpolatlon  and  inverse  FFT  are  Chen  used  Co  compuCe  Che  Cime  domain 
response . 


The  cubic  spline  roucine  used  (CSAKM)  is  a  pare  of  Che  IMSL  Mach  library 
[17] .  Ic  is  based  on  a  meChod  developed  by  Akima  and  is  designed  Co  preserve  che 
shape  of  Che  daca  and  Co  minimize  oscillacions.  The  incerpolacion  is  performed 
on  Che  magnicude  and  phase  arrays,  which  have  been  sorced  by  arranging  Che 
frequencies  in  ascending  order. 


3.6  SMALL  STRUCTURES 

For  small  sCrucCures,  ie.  small  compared  Co  che  wavelengch  of  che  highesc 
frequency  of  che  exciCaCion,  che  firsC  resonance  occurs  above  che  frequency  range 
of  che  excicadon.  Frequency  domain  calculacions  can  be  made  ac  a  relacively  few 
frequencies  and  incerpolaCed  aC  che  incermediaCe  ones  because  chere  are  no  abrupc 
changes  below  che  firsc  resonance. 

Ic  was  observed  on  Figure  3-6  ChaC  a  significanc  pare  of  che  energy  of  che 
EMP  waveforms  lies  on  che  flaC  porcion  of  che  curve,  below  1  MHz.  When 
incerpolacing  Y(w)  Co  perform  che  inverse  cransform,  Af  musc  be  small  enough  co 
cake  chis  porcion  of  che  speccrum  inco  accounc. 
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Figure  3-10  Frequency  domain  response  at  Af  -  0.5  MHz,  obtained  by  interpolating 
the  response  at  varying  Af. 
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Figure  3-11  Solution  obtained  from  standard  application  of  the  FFT  (using 
500  frequencies  up  to  125  MHz)  compared  with  the  solution  (slightly 
shifted)  obtained  by  interpolating  150  frequencies  (up  to  125  MHz. 
variable  Af ) . 
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4.0  CONCLUSIONS 


Ic  was  shown  chat  the  use  of  numerical  eleccromagneclc  Inceracclon 
simulation  codes  for  predicting  the  response  of  complex  systems  to  an  incident 
EMP  field  is  essential  to  the  EMP  analysis  and  hardening  process.  Although  a 
variety  of  codes  exist,  both  in  time  and  frequency  domain,  one  is  often  forced 
to  use  a  time  domain  code  due  to  the  features  or  limitations  inherent  of  each 
code . 


A  number  of  problems  were  solved  with  both  time  domain  and  frequency  domain 
codes.  Ic  was  demonstrated  that  the  use  of  a  frequency  domain  code  can  yield  an 
accurate  solution.  Furthermore,  very  good  agreement  between  different  codes  for 
solving  a  given  problem  was  obtained. 

It  was  also  shown  that  the  old  NATO  definition  of  the  EMP  waveform  (double- 
exponential)  is  not  appropriate  and  that  the  new  definition  (reciprocal)  should 
be  used. 

Techniques  were  developed  to  reduce  the  excessive  computer  time  required 
to  obtain  a  time  domain  solution.  It  was  shown  that  adjusting  the  model  for 
different  frequencies  will  save  about  half  of  the  computer  time.  Ic  was  also 
shown  chat  a  judicious  choice  of  the  frequencies  is  important  to  produce  an 
accurate  solution  with  minimum  computer  time.  A  frequency  domain  interpolation 
technique  was  developed  to  further  reduce  the  computer  time.  All  the  techniques 
developed  can  reduce  the  computer  time  by  a  factor  of  7  to  10. 
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Subroutlna 

Description: 


FFT  (  A.  M  ) 

Radlx-2,  in  place  Fast  Fourier  Transform 
To  perform  Discrete  Fourier  Transform 

This  is  Che  original  algorithm  from  Cooley  &  Tukey 
taken  from  subroutine  HARMM  with  Che  following  changes: 
-  Converted  to  Fortran- 7 7 

•  Use  of  entry  points  instead  of  control  arguments 

•  Scaling  and  pnase  corrected  to  conform  with  the 
definition  of  Che  DFT  algorithm 


Entry  points; 


Arguments : 


FFT  (  A,  M  ) 

FFTS  (  A,  M.  DT  ) 

IFFT  (  A.  M  ) 

IFFTS  (  A.  M.  DF  ) 


Direct  FFT 

Direct  FFT,  with  scaling 
Inverse  FFT 

Inverse  FFT,  with  scaling 


A(*)  Complex  array  of  data,  input  &  output 
M  Log2  of  size  (size  is  2^H) 

DT,  DF  Time  or  frequency  step 


t*-*********** 


Implicit 

Parameter 

Complex 


Rea 


;:j! 


***** 


***** 


Logical 

Entry 

FFT  (  ( 

InvFFc 

-  .False 

Scale 

-  1. 

Goto  10 

Entry 

FFTS  (  , 

InvFFt 

-  .False 

Scale 

-  DT 

Goto  10 

Entry 

IFFT  ( 

InvFFc 

-  .True. 

N 

-  2**M 

Scale 

-  l./N 

Goto  11 

Entry 

IFFTS  ( 

InvFFt 

-  .True. 

Scale 

-  DF 

Incegar*4  (I-M) 

(PI-  3.1415926535) 

A(*),  U,  W.  T 
DT,  DF,  Scale 
InvFFT 


Doing  an  IFFT  ? 
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AFFEMOZX  B 


G^ISELZNES  SUMMARY 


Guldallnes  developed  In  this  report  are  suinmarlzed  herein.  The  primary 
objective  was  to  minimize  the  CPU  time  required  to  run  the  interaction  codes  and 
still  obtain  an  accurate  solution. 


B.l  GUIDELINES  FOR  SOLVINS  TIME  DOMAIN  PROBLEMS 

A  ntimerical  electromagnetic  interaction  code  is  used  to  compute  the 
currents  or  charges  at  multiple  discrete  frequencies.  These  frequencies  need  not 
be  evenly  spaced  and  interpolation  can  be  used  to  obtain  an  evenly  spaced 
sequence  necessary  to  compute  the  inverse  Fourier  transform.  Judicious  choice 
of  frequency  step  and  frequency  range  may  considerably  reduce  the  computer  time. 

a)  Use  a  realistic  model  for  the  excitation  with  a  known  expression  of  its 
Fourier  transform,  such  as  equation  (15). 

b)  Use  longer  segments  and  larger  patches  at  lower  frequencies  to  reduce  the 
problem  size. 

c)  Estimate  the  first  resonance  peak  from  the  dimensions  and  shape  of  the 
object.  Start  with  Af  equal  to  approximately  1/10^^  of  the  estimated  first 
resonance . 

d)  An  incomplete  time  domain  solution,  y(t),  suggests  the  need  for  a  smaller 
Af,  especially  near  the  resonance  peaks. 

e)  The  higher  frequencies  of  Y(b>)  will  contribute  less  and  less  to  the 
solution.  Therefore,  it  is  possible  to  identify  that  trend  in  the  product 
H(u)-E(w)  to  estimate  a  proper  value  for  fa„.  Very  good  results  are 
obtained  when  the  higher  frequencies  contribute  less  chan  1%  of  the  peak 
resonance . 

f)  The  Af  used  to  compute  Che  inverse  transform  (after  interpolation)  should 
be  small  enough  to  cover  Che  lower  frequency  characteristics  (flat  portion 
of  the  spectrum  on  Figure  3-6)  of  the  excitation.  For  the  standard  NATO 
EMF  waveform  definition,  Af  should  be  smaller  chan  1  MH-,;. 


B.2  GUIDELINES  FOR  INTERPOLATING  THE  FREQUENCY  DOMAIN  TRANSFER  FUNCTION 

A  cubic  spline  routine  (CSAKM)  from  the  IMSL  Math  library  is  used  to 
interpolate  the  frequency  domain  solution  computed  at  unevenly  spaced  frequencies 
and  obtain  a  sequence  at  evenly  spaced  frequencies. 


a) 


Convert  the  complex  sequence  into  two  real  sequences  representing  its 
magnitude  and  phase. 


b)  If  necessary,  add  the  DC  value  (£>0).  In  this  case,  Che  magnitude  will  be 
0  and  the  phase  ±90*. 

c)  It  is  necessary  to  unwrap  the  phase  before  Che  interpolation  can  cake 
place . 

d)  Perform  the  interpolation  and  convert  back  to  real  and  imaginary  arrays. 


B.3  GniDEUNES  FOR  PERFORMING  THE  DFT  AND  FFT  ALGORITHMS 

The  application  of  the  DFT  algorithm  requires  chat  Che  N  point  complex 

sequence  be  calculated  at  evenly  spaced  frequencies. 

a)  Multiply  Che  sequence  with  a  window,  such  as  a  cosine* capered  window,  if 
a  significant  discontinuity  is  observed  at  Che  higher  frequencies  because 
of  Che  limited  number  of  frequencies  computed. 

b)  Extend  a  N  point  sequence  to  H  point  by  adding  zeroes  at  the  end.  This 
produces  a  better  interpolation  of  the  solution,  which  is  particularly 
useful  for  plotting,  and  also  allows  Che  use  of  Che  more  efficient  FFT 
algorithm  which  may  have  restriction  on  Che  size  of  che  sequence. 

c)  Mirror  che  M  point  sequence  with  its  conjugate,  thus  forming  a  2M  point 
sequence,  on  which  the  FFT  will  be  performed.  The  real  part  of  result 
will  be  che  solution  (che  imaginary  part  will  be  zero). 

d)  Include  che  effect  of  che  scaling  property  of  the  Fourier  transform  as 
defined  in  equation  (9)  into  che  FFT  algorithm. 

e)  Use  the  equations  (10)  and  (11)  to  compute  At  and  t„,,  from  Af  and  f„„- 
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